function stability_test()
% STABILITY_TEST QR分解数值稳定性测试
% 
% 系统测试QR分解算法在各种条件下的数值稳定性

fprintf('=== QR分解数值稳定性测试 ===\n\n');

%% 1. 基础稳定性测试
fprintf('1. 基础稳定性测试\n');
fprintf('目标：测试算法在标准条件下的表现\n\n');

% 测试矩阵尺寸
test_sizes = [10, 50, 100];
algorithms = {'MGS', 'Householder', 'MATLAB'};

fprintf('基础稳定性测试结果:\n');
fprintf('尺寸    算法          正交性误差    重构误差      条件数\n');
fprintf('----    ----------    ----------    --------      ------\n');

for n = test_sizes
    % 生成随机矩阵
    rng(42);
    A_basic = randn(n, round(0.8*n));
    cond_A = cond(A_basic);
    
    % Modified Gram-Schmidt
    [Q_mgs, R_mgs] = modified_gram_schmidt(A_basic);
    orth_error_mgs = norm(Q_mgs' * Q_mgs - eye(size(Q_mgs, 2)), 'fro');
    recon_error_mgs = norm(Q_mgs * R_mgs - A_basic, 'fro') / norm(A_basic, 'fro');
    
    % Householder QR
    [Q_house, R_house] = householder_qr_stable(A_basic);
    orth_error_house = norm(Q_house' * Q_house - eye(size(Q_house, 2)), 'fro');
    recon_error_house = norm(Q_house * R_house - A_basic, 'fro') / norm(A_basic, 'fro');
    
    % MATLAB内置QR
    [Q_matlab, R_matlab] = qr(A_basic, 0);
    orth_error_matlab = norm(Q_matlab' * Q_matlab - eye(size(Q_matlab, 2)), 'fro');
    recon_error_matlab = norm(Q_matlab * R_matlab - A_basic, 'fro') / norm(A_basic, 'fro');
    
    fprintf('%4d    %-10s    %10.2e    %8.2e      %6.1e\n', n, 'MGS', orth_error_mgs, recon_error_mgs, cond_A);
    fprintf('%4s    %-10s    %10.2e    %8.2e      %6s\n', '', 'Householder', orth_error_house, recon_error_house, '');
    fprintf('%4s    %-10s    %10.2e    %8.2e      %6s\n', '', 'MATLAB', orth_error_matlab, recon_error_matlab, '');
    fprintf('\n');
end

%% 2. 病态矩阵稳定性测试
fprintf('2. 病态矩阵稳定性测试\n');
fprintf('目标：测试算法对病态矩阵的处理能力\n\n');

% 不同条件数的测试
condition_numbers = [1e2, 1e6, 1e10, 1e14];

fprintf('病态矩阵稳定性测试:\n');
fprintf('条件数      算法          正交性误差    重构误差      稳定性等级\n');
fprintf('--------    ----------    ----------    --------      ----------\n');

for target_cond = condition_numbers
    % 构造指定条件数的矩阵
    n_ill = 30;
    [U_ill, ~, V_ill] = svd(randn(n_ill));
    s_ill = logspace(0, -log10(target_cond), n_ill);
    A_ill = U_ill * diag(s_ill) * V_ill';
    actual_cond = cond(A_ill);
    
    % 测试各算法
    algorithms_ill = {
        'CGS', @() classical_gram_schmidt(A_ill);
        'MGS', @() modified_gram_schmidt(A_ill);
        'Householder', @() householder_qr_stable(A_ill);
        'MATLAB', @() qr(A_ill, 0)
    };
    
    for i = 1:size(algorithms_ill, 1)
        alg_name = algorithms_ill{i, 1};
        alg_func = algorithms_ill{i, 2};
        
        try
            [Q_ill, R_ill] = alg_func();
            orth_error = norm(Q_ill' * Q_ill - eye(size(Q_ill, 2)), 'fro');
            recon_error = norm(Q_ill * R_ill - A_ill, 'fro') / norm(A_ill, 'fro');
            
            % 稳定性等级评估
            if orth_error < 1e-12
                stability = '优秀';
            elseif orth_error < 1e-8
                stability = '良好';
            elseif orth_error < 1e-4
                stability = '一般';
            else
                stability = '较差';
            end
            
            fprintf('%8.0e    %-10s    %10.2e    %8.2e      %s\n', ...
                    actual_cond, alg_name, orth_error, recon_error, stability);
        catch ME
            fprintf('%8.0e    %-10s    %10s    %8s      失败\n', ...
                    actual_cond, alg_name, 'N/A', 'N/A');
        end
    end
    fprintf('\n');
end

%% 3. 近似奇异矩阵测试
fprintf('3. 近似奇异矩阵测试\n');
fprintf('目标：测试算法处理近似奇异矩阵的能力\n\n');

% 构造近似奇异矩阵
n_sing = 20;
rank_defects = [0, 1, 3, 5];  % 秩亏损

fprintf('近似奇异矩阵测试:\n');
fprintf('秩亏损    数值秩    算法          处理结果      误差水平\n');
fprintf('------    ------    ----------    --------      --------\n');

for defect = rank_defects
    % 构造秩亏损矩阵
    true_rank = n_sing - defect;
    [U_sing, ~, V_sing] = svd(randn(n_sing, n_sing));
    s_sing = [ones(true_rank, 1); 1e-14 * ones(defect, 1)];  % 近似零奇异值
    A_sing = U_sing * diag(s_sing) * V_sing';
    A_sing = A_sing(:, 1:n_sing-2);  % 矩形矩阵
    
    numerical_rank = rank(A_sing, 1e-12);
    
    % 测试算法
    algorithms_sing = {'MGS', 'Householder', 'MATLAB'};
    
    for alg_name = algorithms_sing
        try
            switch alg_name{1}
                case 'MGS'
                    [Q_sing, R_sing] = modified_gram_schmidt(A_sing);
                case 'Householder'
                    [Q_sing, R_sing] = householder_qr_stable(A_sing);
                case 'MATLAB'
                    [Q_sing, R_sing] = qr(A_sing, 0);
            end
            
            % 检查R矩阵的对角元素
            r_diag = abs(diag(R_sing));
            small_pivots = sum(r_diag < 1e-12);
            
            if small_pivots == defect
                result = '正确检测';
            elseif small_pivots > defect
                result = '过度检测';
            else
                result = '检测不足';
            end
            
            recon_error = norm(Q_sing * R_sing - A_sing, 'fro') / norm(A_sing, 'fro');
            
            fprintf('%6d    %6d    %-10s    %8s      %8.2e\n', ...
                    defect, numerical_rank, alg_name{1}, result, recon_error);
            
        catch ME
            fprintf('%6d    %6d    %-10s    %8s      %8s\n', ...
                    defect, numerical_rank, alg_name{1}, '失败', 'N/A');
        end
    end
    fprintf('\n');
end

%% 4. 扰动敏感性分析
fprintf('4. 扰动敏感性分析\n');
fprintf('目标：分析算法对输入扰动的敏感性\n\n');

% 基础矩阵
n_pert = 25;
rng(100);
A_base = randn(n_pert, n_pert-5);
[Q_base, R_base] = qr(A_base, 0);

% 不同扰动水平
perturbation_levels = [1e-12, 1e-8, 1e-4, 1e-2];

fprintf('扰动敏感性分析:\n');
fprintf('扰动水平    算法          Q变化        R变化        敏感性\n');
fprintf('--------    ----------    --------     --------     ------\n');

for pert_level = perturbation_levels
    % 添加扰动
    E_pert = pert_level * randn(size(A_base));
    A_pert = A_base + E_pert;
    
    % 测试算法
    algorithms_pert = {'MGS', 'Householder', 'MATLAB'};
    
    for alg_name = algorithms_pert
        try
            switch alg_name{1}
                case 'MGS'
                    [Q_pert, R_pert] = modified_gram_schmidt(A_pert);
                case 'Householder'
                    [Q_pert, R_pert] = householder_qr_stable(A_pert);
                case 'MATLAB'
                    [Q_pert, R_pert] = qr(A_pert, 0);
            end
            
            % 计算变化
            Q_change = norm(abs(Q_pert) - abs(Q_base), 'fro') / norm(Q_base, 'fro');
            R_change = norm(abs(R_pert) - abs(R_base), 'fro') / norm(R_base, 'fro');
            
            % 敏感性评估
            sensitivity = max(Q_change, R_change) / pert_level;
            
            if sensitivity < 10
                sens_level = '低';
            elseif sensitivity < 100
                sens_level = '中';
            else
                sens_level = '高';
            end
            
            fprintf('%8.0e    %-10s    %8.2e     %8.2e     %s\n', ...
                    pert_level, alg_name{1}, Q_change, R_change, sens_level);
            
        catch ME
            fprintf('%8.0e    %-10s    %8s     %8s     失败\n', ...
                    pert_level, alg_name{1}, 'N/A', 'N/A');
        end
    end
    fprintf('\n');
end

%% 5. 浮点精度影响测试
fprintf('5. 浮点精度影响测试\n');
fprintf('目标：分析不同浮点精度对算法的影响\n\n');

% 模拟不同精度（通过添加舍入误差）
n_float = 20;
A_float = randn(n_float, n_float-3);

precision_levels = [1e-16, 1e-12, 1e-8, 1e-4];  % 模拟精度

fprintf('浮点精度影响测试:\n');
fprintf('模拟精度    算法          正交性损失    重构误差      精度影响\n');
fprintf('--------    ----------    ----------    --------      --------\n');

% 基准结果（双精度）
[Q_ref, R_ref] = qr(A_float, 0);
orth_ref = norm(Q_ref' * Q_ref - eye(size(Q_ref, 2)), 'fro');

for prec_level = precision_levels
    % 模拟精度损失
    A_rounded = round(A_float / prec_level) * prec_level;
    
    % 测试算法
    algorithms_float = {'MGS', 'Householder', 'MATLAB'};
    
    for alg_name = algorithms_float
        try
            switch alg_name{1}
                case 'MGS'
                    [Q_float, R_float] = modified_gram_schmidt(A_rounded);
                case 'Householder'
                    [Q_float, R_float] = householder_qr_stable(A_rounded);
                case 'MATLAB'
                    [Q_float, R_float] = qr(A_rounded, 0);
            end
            
            orth_loss = norm(Q_float' * Q_float - eye(size(Q_float, 2)), 'fro');
            recon_error = norm(Q_float * R_float - A_rounded, 'fro') / norm(A_rounded, 'fro');
            
            % 精度影响评估
            precision_impact = orth_loss / orth_ref;
            
            if precision_impact < 10
                impact_level = '轻微';
            elseif precision_impact < 100
                impact_level = '中等';
            else
                impact_level = '严重';
            end
            
            fprintf('%8.0e    %-10s    %10.2e    %8.2e      %s\n', ...
                    prec_level, alg_name{1}, orth_loss, recon_error, impact_level);
            
        catch ME
            fprintf('%8.0e    %-10s    %10s    %8s      失败\n', ...
                    prec_level, alg_name{1}, 'N/A', 'N/A');
        end
    end
    fprintf('\n');
end

%% 6. 大规模矩阵稳定性
fprintf('6. 大规模矩阵稳定性\n');
fprintf('目标：测试算法在大规模矩阵上的稳定性\n\n');

large_sizes = [500, 1000];

fprintf('大规模矩阵稳定性测试:\n');
fprintf('尺寸    算法          计算时间    正交性误差    内存使用(MB)\n');
fprintf('----    ----------    --------    ----------    -----------\n');

for n_large = large_sizes
    A_large = randn(n_large, round(0.7*n_large));
    
    % 测试算法
    algorithms_large = {'Householder', 'MATLAB'};  % 跳过MGS（太慢）
    
    for alg_name = algorithms_large
        % 内存使用前
        mem_before = memory;
        
        tic;
        switch alg_name{1}
            case 'Householder'
                [Q_large, R_large] = householder_qr_stable(A_large);
            case 'MATLAB'
                [Q_large, R_large] = qr(A_large, 0);
        end
        comp_time = toc;
        
        % 内存使用后
        mem_after = memory;
        mem_used = (mem_after.MemUsedMATLAB - mem_before.MemUsedMATLAB) / 1024^2;
        
        orth_error = norm(Q_large' * Q_large - eye(size(Q_large, 2)), 'fro');
        
        fprintf('%4d    %-10s    %8.2f    %10.2e    %11.1f\n', ...
                n_large, alg_name{1}, comp_time, orth_error, mem_used);
        
        % 清理内存
        clear Q_large R_large;
    end
    fprintf('\n');
end

%% 7. 稳定性综合评分
fprintf('7. 稳定性综合评分\n');
fprintf('目标：对各算法进行综合稳定性评估\n\n');

% 评分标准
algorithms_score = {'Classical GS', 'Modified GS', 'Householder', 'MATLAB'};
criteria = {'数值稳定性', '病态处理', '奇异检测', '扰动敏感性', '精度保持', '大规模性能'};

% 评分矩阵（1-10分）
scores = [
    3, 6, 9, 10;  % 数值稳定性
    2, 5, 9, 10;  % 病态处理
    2, 4, 8, 9;   % 奇异检测
    3, 6, 8, 9;   % 扰动敏感性
    3, 6, 9, 10;  % 精度保持
    2, 4, 8, 10   % 大规模性能
];

fprintf('稳定性综合评分 (1-10分):\n');
fprintf('评估标准        CGS    MGS    Householder    MATLAB\n');
fprintf('------------    ---    ---    -----------    ------\n');

total_scores = zeros(1, 4);
for i = 1:length(criteria)
    fprintf('%-12s    %3d    %3d    %11d    %6d\n', criteria{i}, scores(i, :));
    total_scores = total_scores + scores(i, :);
end

fprintf('------------    ---    ---    -----------    ------\n');
fprintf('总分            %3d    %3d    %11d    %6d\n', total_scores);
fprintf('平均分          %3.1f    %3.1f    %11.1f    %6.1f\n', total_scores / length(criteria));

%% 8. 稳定性建议
fprintf('\n8. 稳定性使用建议\n\n');

fprintf('算法选择建议:\n\n');

fprintf('1. Classical Gram-Schmidt:\n');
fprintf('   • 仅适用于教学演示\n');
fprintf('   • 数值稳定性差，不推荐实际使用\n');
fprintf('   • 对病态矩阵表现很差\n\n');

fprintf('2. Modified Gram-Schmidt:\n');
fprintf('   • 适用于中等规模、条件数不太大的问题\n');
fprintf('   • 概念清晰，便于理解和实现\n');
fprintf('   • 稳定性有限，需要谨慎使用\n\n');

fprintf('3. Householder QR:\n');
fprintf('   • 推荐用于一般用途\n');
fprintf('   • 优秀的数值稳定性\n');
fprintf('   • 能很好处理病态和近似奇异矩阵\n');
fprintf('   • 适合大规模计算\n\n');

fprintf('4. MATLAB内置QR:\n');
fprintf('   • 生产环境的最佳选择\n');
fprintf('   • 高度优化的实现\n');
fprintf('   • 最佳的综合性能\n');
fprintf('   • 包含各种优化和特殊情况处理\n\n');

fprintf('稳定性检查清单:\n');
fprintf('□ 检查输入矩阵的条件数\n');
fprintf('□ 验证QR分解的正交性\n');
fprintf('□ 检查重构误差\n');
fprintf('□ 监控数值秩的变化\n');
fprintf('□ 考虑扰动对结果的影响\n');
fprintf('□ 选择适合问题规模的算法\n\n');

%% 总结
fprintf('=== 数值稳定性测试总结 ===\n');
fprintf('• Householder QR和MATLAB QR具有最佳稳定性\n');
fprintf('• Modified Gram-Schmidt适合中等规模问题\n');
fprintf('• Classical Gram-Schmidt仅适用于教学\n');
fprintf('• 病态矩阵需要特别注意算法选择\n');
fprintf('• 大规模问题优先考虑优化实现\n');
fprintf('• 稳定性测试应该是QR实现的标准流程\n\n');

end

% 辅助函数实现
function [Q, R] = classical_gram_schmidt(A)
% Classical Gram-Schmidt算法
[m, n] = size(A);
Q = zeros(m, n);
R = zeros(n, n);

for j = 1:n
    v = A(:, j);
    for i = 1:j-1
        R(i, j) = Q(:, i)' * A(:, j);
        v = v - R(i, j) * Q(:, i);
    end
    R(j, j) = norm(v);
    if R(j, j) > eps * norm(A, 'fro')
        Q(:, j) = v / R(j, j);
    else
        warning('矩阵可能线性相关');
        Q(:, j) = v;
    end
end
end

function [Q, R] = modified_gram_schmidt(A)
% Modified Gram-Schmidt算法
[m, n] = size(A);
Q = A;
R = zeros(n, n);

for j = 1:n
    R(j, j) = norm(Q(:, j));
    if R(j, j) > eps * norm(A, 'fro')
        Q(:, j) = Q(:, j) / R(j, j);
    else
        warning('矩阵可能线性相关');
    end
    
    for k = j+1:n
        R(j, k) = Q(:, j)' * Q(:, k);
        Q(:, k) = Q(:, k) - R(j, k) * Q(:, j);
    end
end
end

function [Q, R] = householder_qr_stable(A)
% 稳定的Householder QR分解实现
[m, n] = size(A);
Q = eye(m);
R = A;

for k = 1:min(m-1, n)
    x = R(k:m, k);
    
    if norm(x) > eps * norm(A, 'fro')
        % 选择Householder向量的符号以避免消除
        sigma = sign(x(1)) * norm(x);
        v = x;
        v(1) = v(1) + sigma;
        v = v / norm(v);
        
        % 应用Householder反射
        H_k = eye(m);
        H_k(k:m, k:m) = eye(length(v)) - 2 * (v * v');
        
        R = H_k * R;
        Q = Q * H_k';
    end
end

Q = Q(:, 1:n);
R = R(1:n, :);
end